#include <vector>


struct Pos
	{
	int x, y;
	};

template<class CONT, class EL>
void remove(CONT & c, const EL & e)
	{
	for (typename CONT::iterator i = c.begin(); i!= c.end(); i++)
		if (*i == e)
			{
			*i = c.back();
			c.pop_back();
			return;
			}
	}

struct Z
	{
	unsigned z1, z2, z3, z4;
	};

// S1, S2, S3, and M are all constants, and z is part of the  
// private per-thread generator state.  
unsigned tausStep(unsigned &z, int S1, int S2, int S3, unsigned M)  
	{  
//	const unsigned b=(((z << S1) ^ z) >> S2);  
	return z = (((z & M) << S3) ^ (((z << S1) ^ z) >> S2));  
	}

unsigned lcgStep(unsigned & z, int A, unsigned C)
	{
	z = A * z + C;
	}

unsigned hybridTausU(Z & z)  
	{  
	// combined period is lcm(p1,p2,p3,p4)~ 2^121  
	return   
		tausStep(z.z1, 13, 19, 12, 4294967294ul) ^  // p1=2^31-1  
		tausStep(z.z2, 2, 25, 4, 4294967288ul) ^    // p2=2^30-1  
		tausStep(z.z3, 3, 11, 17, 4294967280ul) ^   // p3=2^28-1  
		lcgStep(z.z4, 1664525, 1013904223ul)        // p4=2^32  
		;  
	}  

const int nPop = 256;
const int nInd = 256;
const int maxX = 1024;
const int maxY = 1024;
const int cellSz = 64;
const int cacheMX = maxX / cellSz;
const int cacheMY = maxY / cellSz;
const int vis = 10;
const int maxSqDist = vis * vis;

std::size_t cachePos(std::size_t x)
	{
	return x / cellSz;
	}

typedef std::vector<std::size_t> cell_t;

int main()
	{
	Z z;
	z.z1 = 5093;
	z.z2 = 384938;
	z.z3 = 93010;
	z.z4 = 34384;

	std::vector<std::vector<Pos> > pos;
	std::vector<std::vector<std::vector<cell_t> > > cache;
	std::vector<std::vector<int> > coll;

	pos.resize(nPop);
	cache.resize(nPop);
	coll.resize(nPop);

	for (std::size_t p=0; p<pos.size(); p++)
		{
		pos[p].resize(nInd);
		coll[p].resize(nInd);
		cache[p].resize(cacheMX);
		for (std::size_t  x=0; x<cacheMX; x++)
			cache[p][x].resize(cacheMY);

		for (std::size_t i=0; i<pos[p].size(); i++)
			{
			pos[p][i].x = hybridTausU(z) % maxX;
			pos[p][i].y = hybridTausU(z) % maxY;
			cache[p][cachePos(pos[p][i].x)][cachePos(pos[p][i].y)].push_back(i);
			}
		}

	for (std::size_t s=0; s<1000; s++)
		{
		for (std::size_t p=0; p<pos.size(); p++)
			{
			for (std::size_t i=0; i<pos[p].size(); i++)
				{
				const int cpX = cachePos(pos[p][i].x);
				const int cpY = cachePos(pos[p][i].y);

				const int newX = (pos[p][i].x + hybridTausU(z)%5 - 2 + maxX) % maxX;
				const int newY = (pos[p][i].y + hybridTausU(z)%5 - 2 + maxY) % maxY;

				const int cpnX = cachePos(newX);
				const int cpnY = cachePos(newY);

				if (cpnX != cpX || cpnY != cpY)
					{
					remove(cache[p][cpX][cpY], i);
					cache[p][cpnX][cpnY].push_back(i);
					}
				
				pos[p][i].x = newX;
				pos[p][i].y = newY;
				}
			for (std::size_t i=0; i<pos[p].size(); i++)
				{
				const int minCpX = cachePos(pos[p][i].x - vis);
				const int minCpY = cachePos(pos[p][i].y - vis);
				const int maxCpX = cachePos(pos[p][i].x + vis);
				const int maxCpY = cachePos(pos[p][i].y + vis);

				for (int cpx=minCpX; cpx<=maxCpX; cpx++)
					{
					const int offs_x = cpx<0 ? -maxX : 
						(cpx>=cacheMX ? maxX : 0);
					const int rx = (cpx+cacheMX)%cacheMX;

					for (int cpy=minCpY; cpy<=maxCpY; cpy++)
						{
						const int offs_y = cpy<0 ? -maxY : 
							(cpy>=cacheMY ? maxY : 0);
						const cell_t & local = cache[p][rx][(cpy+cacheMY)%cacheMY];
						for (cell_t::const_iterator l=local.begin(); 
							l!=local.end(); l++)
							{
							const int dx = pos[p][i].x - (pos[p][*l].x + offs_x);
							const int dy = pos[p][i].y - (pos[p][*l].y + offs_y);
							const int sqDist = dx*dx + dy*dy;
							if (sqDist < maxSqDist)
								coll[p][i]++;
							}
						}
					}
				}
			}
		}
	}
